Corticothalamic projections control synchronization in locally coupled bistable thalamic oscillators 
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Thalamic circuits are able to generate state-dependent oscillations of different frequencies and degrees of 
synchronization. However, only little is known how synchronous oscillations, like spindle oscillations in the 
thalamus, are organized in the intact brain. Experimental findings suggest that the simultaneous occurrence of 
spindle oscillations over widespread territories of the thalamus is due to the corticothalamic projections, as the 
synchrony is lost in the decorticated thalamus. Here we study the influence of corticothalamic projections on the 
synchrony in a thalamic network, and uncover the underlying control mechanism, leading to a control method 
which is applicable in wide range of stochastic driven excitable units. 
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gins and mechanisms of sleep and other processes in the mam- 
malian brain, and - in the long run - to control them. 

The paper is organized as follows: We first describe our 
computational model, which is a simplified phenomenolog- 
ical model of the thalamocortical oscillator described de- 
Then, in a first startup, we study burst- 



Coupled oscillators are abundant in physics UJ, |2|], chem- 
istry yj and biology S 01 ■ Whenever large numbers of 
coupled oscillators are considered-the collective behavior of 
the ensemble is of great interest yl 01- A widespread phe- 
nomenon in populations of periodic, noisy and chaotic oscil- 
lators (or maps) is the appearence of synchronous collective 
oscillations, studied theoretically [7] as well as experimen- 
tally Jj, |8[]. In neural systems the phenomenon of spike burst 
activity is widespread and of great importance |9lll0|]. This 
activity is characterized by a recurrent transition between a 
resting state and a firing state with a burst of multiple spikes. 
Bursting is a multiple time scale effect and appears due to a 
slow process which modulates a fast subsystem lr9l [Tlll . 

Neural cells in the thalamus are known to exhibit spike 
burst activity during periods of drowsiness, inattentiveness 
and sleep 11211 . Spindle oscillations, a hallmark of early 
sleep stages, observed in the electroencephalograhic record- 
ings of sleeping humans as oscillations in a 12-15 Hz fre- 
quency range, are considered to be the result of synchronized 
spike burst activity of millions of neurons in the thalamus 
I loL 131. Despite experimental findings in thalamic slices, 
where the spindle oscillations propagate like a traveling wave 
through the network, in the intact brain spindle oscillations 
occur almost synchronously over widespread territories of the 
thalamus 113, 11411 . These contrasting results between in vitro 
and in vivo experiments occur due to the absence of corti- 
cothalamic projections (synaptic connections from the cortex 
to the thalamus) in a thalamic slice as, after ablation of the 
cortex, synchrony is lost in the thalamus lll3lll4ill5ll . 

In a coupled system of oscillators, oscillations become 
more and more coherent when the coupling strength is in- 
creased ||2|, lZ|]. Here we explore the possibility to induce 
the transition from the decoherent state to the coherent state 
- in which bursts occur synchronously - by external signals. 
For this reason, we investigate the mechanism by which cor- 
ticothalamic projections control coherence of thalamic spin- 
dle oscillations. The knowledge of control mechanisms for 
synchrony in neural systems may help to understand the ori- 



tailed in [12, 

synchronization analytically for the case of two coupled neu- 
rons. Finally, we study the collective behavior of thalamic 
oscillators numerically, and verify the experimental observa- 
tion that corticothalamic projections control synchrony in the 
thalamus by using a hybrid network, consisting of a computa- 
tional model of the thalamus and real human slow wave EEG 
data as the control signal. 

The dynamics of the system described in 11611 is mainly de- 
termined by a slow change between a bistable state where 
a stable fixed point and a stable limit cycle coexist, and 
a monostable state where only a stable fixed point exists. 
This behavior can be described by a slowly modulated over- 
damped movement of a particle in a rotating symmetric time- 
dependent double well potential (Fig.[TJ 
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where r — x 2 + y 2 , and (pit) = tot captures the oscil- 
lation in phase (mod 2tt), see Ref. 1 17] for a related model. 
Depending on the slow variable a(t), the system possess ei- 
ther one stable fixed point, or it shows bistability with a sta- 
ble fixed point and a stable limit-cycle. In this paper, we re- 
strict ourselves to very slow changes in ait). In that case, 
the relaxation time within the wells is short compared to ait), 
so the quasi-adiabatic approximation is applicable. We get 
a Langevin equation for the movement in r direction for the 
single oscillator: 



dU(r, a) 
dr 



Fit) 



u = 0.01(r(l -u)- 0.14u), (2) 
where ait) = a(l — u), F(t) is a common control signal. 
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FIG. 1: Left: Depending on a, the system possesses either one sta- 
ble fixed point or it shows bistability with a stable fixed point and a 
stable limit-cycle. Right: Bifurcation diagram of eq. ([2} with a(t) 
considered as a bifurcation parameter. Solid lines indicate a stable 
fixed point, open circles an unstable limit cycle and filled circles a 
stable limit cycle. Depending on a, the system possesses either one 
stable fixed point, or it shows bistability with a stable fixed point and 
a stable limit cycle. The limit-cycle dissapears by an inverse sniper 
(saddle node into periodic orbit) bifurcation. 

Following [9], as a first step we study synchrony in two 
bistable elements coupled linearly via the r component: 

fi = -n + a(l - u t )r 3 - r\ + If xt + etj 

v-i = - u i) - 0.14uj), (3) 

where e > is the coupling strength, a = 3.1 is chosen such 
that U (r, a) is a double well potential for u < 0.35, and 
i = 1,2, j = 2, 1 respectively are the indices. If xt is the 
external input, which is a Poisson distributed shotnoise with 
medium rate of 1 /100ms, refractory period of 30ms and pulse 
duration of 2ms. We emphasize that the If xt are stochastically 
independent, so {/f xt , IJ xt ) = <5;.j. The synchronized state is 
then given by 

f = —r + a(l — u)r 3 — r 5 + er 

u = fJ.{r(l — u) — bu), (4) 

as in this case |ri — ra| and — it%\ vanish for t — > 00. 
Now we consider the temporal evolution of small disturbances 
of the synchronous manifold. For this manner we trans- 
form to r± = r\ — T2 and u± = u\ — 112. If the sys- 
tem is close to the synchronous state \r±\, \u±\ <C 1, the 
following approximations hold: r™ — r% = nr n ~ 1 r± and 
r"ui — r% U2 = r n u± + unr n ~ 1 r±. So the movement away 
from the synchronous state can be described by: 

r± = — rj_+ 3ar 2 r±(l — u) — aT 3 u± — 5r i r± + Ij**— er± 
U± = /i((l — u)r± — ru± — bu±), (5) 

The stability of the synchronized state is determined by the 
solution of eqns. ( 14151 . To estimate the coupling strength at 
which burst synchronization occurs, we consider the quasi- 
adiabatic limit fi — > in eqns. i4\5i . what represents the 



r— subsystem: 

r = — r + a(l - u)r 3 - r 5 + e(r + Ar) (6) 
r_L = r ± 7(r, u) + 1^ - ru ± , (7) 

where 7(5", u) = (—1 + 3ar 2 (l — u) — 5r 4 — e), and Ar is 
the deviation from the synchronous state. The synchronized 
state of (0 is only stable, if the Lyapunov exponent of (0 
is negative. Further, in the quasi-adiabatic approximation we 
have du/dt = du±/dt = so that u and u± can be treated 
as parameters in eqns. ||6JQ. In ©, we have to take care of 
the input I cxt which disturbs the synchronous manifold; we 
take respect of this effect by the term A7\ In order to be able 
to estimate the Lyapunov exponent, we need to get r(u). For 
this purpose we use the fact, that due to the geometry of the 
potential U(r, t) in fl}, the only area where small disturbances 
r± <C 1 can grow in time, is the surrounding of the local 
maximum of U(r, t). So we get the following conditions for 
r(u): 

dU/dr = and d 2 U/dr 2 < (8) 

As a condition for the stability of the synchronous manifold 
we have j(r, u) < 0. By setting u = 0.12 in (13), we get r(e) 
for the case of a maximal sensitivity of the synchronous mani- 
fold to disturbances, as in this case the system has its minimal 
threshold. Although the linearization is only valid for r± <C 1 
we will consider the case where one neuron is in the resting 
state and the other neuron is in the excited state, which leads to 
values of r± > 1, in this case our linearization only delivers a 
lower bound estimation for e. Here the approximation is sup- 
ported by the observation that for increasing e, first the jumps 
between the two wells synchronize and then the intra-well os- 
cillations get more and more coherent. Solving j(r, e) > 
for A r = 1.1 and Ar = 0.5 numerically, 7 changes its sign 
from negative to positive at e ~ 0.16 for Ar = 1.1 and at 
e ~ 0.25 for Ar — 0.5. If we compare Fig. |2]we see that our 
estimation works also quite well for the case discussed above. 
Due to this estimation of the coupling strength for the tran- 
sition to the synchronized state, the variations of r± should 
decrease drastically at e — > 0.16 as here the transitions be- 
tween the two wells get synchronized and should show only 
little changes for e > 0.25. The numerical simulations in Fig. 
|2]reproduce the results of our estimation quite well. The esti- 
mation above only gives the values of e where r synchronizes, 
but however as u obeys a linear differential equation, which 
gets activated by r, we expect u to synchronize at the same 
value of e as r. The numerical simulation of eqns. ([3]l for 
different values of e in Fig. |2]verify this argumentation. 

Now we model the thalamic network by a 2-dimensional 
50 x 50 square lattice with a nearest neighbor coupling and 
periodic boundary conditions. The network of N neurons will 
be described by the following equation 
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FIG. 2: a): 7 for the case that the synchronous manifold from 10 gets 
disturbed, by an input spike with amplitude 0.5 in the surrounding of 
the local maximum of U (r, a), see {8]l and Fig. [2 b) 7 for the case 
that the two neurons are in different wells of the potential, c) and 
d): Variations of r± and u± as a function of e. As the interwell 
jumps synchronize at e ~ 0.16, the variations of r± (c) and tii (d) 
rapidly decrease at this point. For e > 0.26 the variations show no 
more significant changes as here the whole synchronous manifold is 
stable. 



where the uncoupled dynamics of the (ij)-th node obeys 



) given by (0 and the N x N matrix G determines 



the coupling between the neurons. For a nearest neighbor cou- 
pling its components are given by Gij — 5i.j+i + with 
periodic boundaries = N. F cxt (i) = n£(t) is the common 
external forcing which consists of real human slow wave sleep 
EEG-data £(t). Almost any coupling scheme can be cast into 
the form of eq. (O by choosing the right G matrix 12011 . here 
we use a linear nearest neighbor coupling with radius 1 with- 
out self-loops. In a first instance, we study the system with- 
out external control, such F(t) = 0. In dependence of the 
coupling strength e we observe four different kinds of collec- 
tive phenomena. For low values of e < 0.029 the oscillators 
are desynchronized. Due to the nearest neighbor coupling for 
0.029 < e < 0.057 the bursts propagate like traveling waves 
through the network. For e > 0.057, burst-synchronization 
occurs (see Fig. 01. For even larger e, a trivial homogeneous 
state (not shown) occurs. In this work we are mainly inter- 
ested in the transition of the traveling bursts to synchronous 
network bursts, and the possibility to induce this transition by 
an external control signal. As mentioned above, it was ob- 
served in thalamic slices that spindle oscillations propagate in 
a way similar to traveling waves in the absence of corticotha- 
lamic projections II 1 3L 1 1 411 . In II 1011 it was examined whether 
a comparable temporal grouping of spindle activity, coincid- 
ing with cortical slow wave oscillations can be found during 
human slow wave sleep. The results clearly show that also in 
the human sleep the spindles become synchronized by corti- 
cothalamic projections. This phenomenon was part of several 
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FIG. 3: Depending on e the network shows either asynchronous (e < 
0.029) or synchronous oscillations (e > 0.057). Between these two 
states traveling waves occur. 



numerical and experimental investigations II15LI21H . 

At this point we elucidate for a qualitative understanding 
of the network behaviour. We assume the slow process u^ (t) 
to be a relaxation oscillator with relaxation time t, this re- 
laxation process gets activated by the variable ry (t), in turn 
(t) inhibits excitation of during the relaxation time t. 
So (t) can be interpreted as self-inhibitory process which 
leads to a dead time, in which the single oscillator is not ex- 
citable. So the maximum phase difference between any ar- 
bitrary chosen oscillators is at most of the magnitude of r. 
Another way to block excitation of all is a common suffi- 
ciently strong hyperpolarization of all oscillators, if the dura- 
tion of this hyperpolarization exceeds t, all Uij (t) will decay 
back to their excitable ground state. That means they all have 
the same phase. If we further assume the input to be a Poisson 
distributed spike train with medium rate A, and further assume 
that in the resting state u^ (t) gets activated by a single spike, 
then the mean waiting time until Uy (t) gets activated is given 
by 1/A. So for A -1 <C r the Uij(t) can be synchronized 
by a sufficiently long hyperpolarization. This effect is often 
called phase resetting. Of course the phase of the u^ (t) can 
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perimental results in [10]. Fig. [4] shows that the correlation 
between driving EEG signal and meanfield activity behaves 
similar to the corresponding experimental result of Ref. 1 1 Ofl - 

To conclude, we have studied the synchronizing influence 
of corticothalamic projections in a thalamic network by a the- 
oretical model and computationally reproduced the experi- 
ment. As the model equations (apart from parameter choices) 
do not rely specifically on a neural substrate, we anticipate the 
qualitative behavior to be generic also for other stochastically 
driven systems composed of excitable units. 

This research has been supported by the Deutsche 
Forschungsgemeinschaft (DFG) within SFB 654 "Plasticity 
and Sleep". The authors thank Jan Born and Lisa Marshall 
for intensive and fruitful discussions. 



FIG. 4: Correlation of meanfield network activity with driving input 
signal depending on the driving parameter k for e = 1 /25. a) and c): 
Numerical simulation of the network driven by EEG spindle sleep 
data as used in 1 10]. For k — 0.0 correlations vanish, b) and d): The 
corresponding experimentally obtained correlations between spindle 
rms which corresponds to the of the nj scaled by a factor of 30 and 
EEG activity. 

diffuse by time because of the highly stochastic input. But 
if the phase resetting is done repetitively, then the variance 
of the phases of can be confined. So we post that cor- 
ticothalamic projections control coherence in thalamic slices 
by an open loop repetitive phase resetting. The effect of phase 
resetting to synchronize stochastically excited relaxation os- 
cillators gets enforced, if the hyperpolarization is followed by 
a depolarization, as this increases the probability of excita- 
tion (see Fig [4] a) an d c ))- in this sense, cortical slow waves, 
which consist of a hyperpolarizing half-wave followed by de- 
polarizing half-wave, form an optimized signal to synchronize 
stochastic relaxation oscillators. 

Finally, we verify our approach computationally in compar- 
ison with experimental data. Consequently, we use a hybrid 
network consisting of a computational model of the thalamic 
slice, which gets real slow wave sleep EEG-data as a com- 
mon input, to investigate the control of coherence in the tha- 
lamic network. While one might argue that the architecture 
of this model does not reflect nature in detail, as the thalam- 
ocortical projections are completely neglected, one reason for 
this approximation is the fact that the corticothalamic projec- 
tions outnumber the thalamocortical projections by a number 
of magnitude 11311 : second it was shown in cortical slices that 
slow waves also occur in the isolated cortex i.e. without tha- 
lamic input @]. From both reasons we assume that synchro- 
nization in the thalamus mainly is controlled through open- 
loop control by the cortex. The results obtained by our hybrid 
network provide a minimal model comparable with the ex- 
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